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I.  INTRODUCTION 


The  concept  of  distributing  the  functionality  of  larger  satellites  among  smaller, 
cooperative  satellites  has  been  seriously  considered  for  assorted  space  missions  to 
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accomplish  goals  that  are  not  possible  or  very  difficult  to  do  with  a  single  satellite. 

The  current  trend  in  satellite  design  is  toward  using  smaller,  more  capable  satellites  in 
cooperative  formations  or  distributed  arrays.4'5 

An  extremely  critical  issue  in  architecting  satellite  swarms  is  the  formation 
configuration.  Arguably,  the  main  feasibility  criterion  in  the  architecture  of  satellite 
swarms  is  the  design  of  globally-minimum-fuel  configurations.  This  simply  follows  from 
the  fact  that  fuel  for  orbiting  spacecraft  comes  at  a  very  high  premium  and  could 
significantly  offset  any  other  advantage  held  by  a  swarm  configuration.  Thus,  there  is  a 
need  to  determine  zero-propellant  fonnation  configurations  (if  they  exist)  and  methods 
for  controlling  the  formation  with  little  or  no  propellant.  It  is  well  known  that  a  family  of 
zero-propellant  circular  and  elliptic  formations  exists  when  the  spacecraft  are  subject 
only  to  an  inverse-square  gravity  field.  '  However,  these  formations  tear  apart  in  the 
presence  of  “disturbing”  effects  such  as  J2.  Thus,  a  search  for  invariant  relative  orbits  or 
formations  (if  they  exist)  goes  on.  Another  effect  that  must  be  accounted  for  is  a  non¬ 
zero  eccentricity  of  the  reference  orbit.9  Unlike  the  J2  disturbance,  an  error  in  eccentricity 
can  be  controlled  by  a  one-time  expenditure  of  propellant.  This  is  based  on  the 
observation  that  the  J2  disturbance  is  an  error  in  the  dynamical  model  whereas  the  error  in 
eccentricity  is  one  of  initial  conditions.  The  error  in  eccentricity  does  not  fully  address 
the  problem  since  in  many  applications  it  is  desirable  to  have  formations  for  every 
eccentricity  (and  not  just  small  eccentricities).  Hence,  the  “real”  problem  is  to  find 
formations  in  the  presence  of  the  totality  of  (modelable)  deterministic  forces.  If  such 
formations  do  not  exist  naturally,  then  it  is  imperative  that  the  minimum-fuel  formation 
configuration  be  detennined.  The  total  fuel  consumption  of  a  formation  over  its  lifetime 
then  establishes  the  practical  feasibility  of  a  swarm  architecture.  There  are  some 
opinions10  suggesting  that  minimum-fuel  configurations  exist  and  that  they  do  not  offset 
the  advantages  offered  by  other  perfonnance  metrics.  However,  there  appears  to  be  no 
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general-purpose  method  published  in  the  open  literature  on  how  to  find  these  minimum- 
fuel  configurations  subject  to  arbitrary  forces. 

This  thesis  presents  a  general  problem  algorithm  for  finding  optimal  formations 
and  shows  that  a  very  natural  setting  for  solving  this  problem  is  (nonlinear)  optimal 
control  theory.  In  particular,  these  ideas  are  formulated  by  using  elements  from  optimal 
periodic  control  theory  with  partially  periodic  states .n  The  Legendre  pseudospectral 
technique  is  applied  to  numerically  solve  this  problem  (see  Reference  12  and  the 
references  contained  therein).  Because  of  the  Covector  Mapping  Theorem,  this 
technique  is  neither  a  direct  nor  an  indirect  method.  Rather,  it  provides  all  the  ease  of  a 
direct  method  while  providing  the  accuracy  of  an  indirect  method.  This  method  is 
implemented  using  the  general-purpose  software  package,  DIDO,14  which  has  been  used 
extensively  over  the  past  few  years  to  solve  a  myriad  of  complex  optimal  control 
problems. 

The  applicability  of  optimal  control  theory  to  satellite  fonnation  was  not  known  at 
the  outset  of  the  thesis  work  presented  here.  For  this  reason,  the  research  style  chosen 
was  a  building-block  approach.  First,  optimal  control  theory  had  to  be  demonstrated  as 
an  appropriate  framework  for  a  simplified  model  and  then  complexity  would  be  added 
incrementally. 

A  model  simply  captures  the  essential  aspects  from  a  certain  point  of  view  and 
simplifies  or  ignores  the  rest.  The  so-called  Hill-Clohessy-Wiltshire15  (C-W)  equations 
were  chosen  as  the  first  model  specifically  because  the  solutions  were  known.  This 
allowed  a  validation  of  the  method  before  embarking  upon  more  complicated  models. 
The  second  model  was  chosen  to  address  the  eccentric  reference  orbit.  To  do  this,  a 
nonlinear  version  of  the  linearized  equations 16,1 7,1 8  was  used  for  the  dynamical  model 
without  the  presumption  of  a  solution. 
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II.  GENERAL  FRAMEWORK 


A.  PROBLEM  DEFINITION19 

The  notion  of  a  swarm  can  be  defined  by  stating  that  a  group  of  satellites  are  said 
to  be  in  formation  if  a  given  configuration  metric  c({x'})  is  bounded, 

c,  <c({x'})<c„  (2.1) 

where  x'  denotes  the  state  of  the  ih  spacecraft  at  time  x.  This  state  can  be  the  usual 

position-velocity  state  or  any  other  set  (e.g.  orbital  elements). 

The  dynamics  of  a  swarm  of  Ns  spacecraft  are  given  in  some  coordinate  system 

by, 

x'=f(x\u,x;p)  i  =  l...Ns  (2.2) 

If  the  distance  between  any  two  satellites,  d(x‘  ,xJ)  ,  is  chosen  as  the  metric,  then 
the  swarm  is  said  to  be  in  formation  if 

df  <d{xl,xJ)<df  Vx  (2.3) 

where  di  and  du  define  the  smallest  and  largest  allowable  separation  distances, 
respectively.  Instead  of  choosing  separation  distances  between  every  spacecraft  pair,  it  is 
sometimes  simpler  to  choose  a  separation  distance  between  a  spacecraft  and  a  reference 
spacecraft.  In  this  case,  equation  (2.3)  can  be  replaced  by 

d{  <  d(y,xJ)<d]u  Vx  (2.4) 

where  y  is  the  state  of  the  reference  spacecraft.  From  these  fundamentals,  a  family  of 
formations  can  be  defined  as  follows.  If  d(y,xJ)  is  a  constant  for  all  j,  then  the 
formation  is  a  circular  formation, 

d(y,x])  =  d{  =  dJu  Vx  (2.5) 

with  the  spacecraft  at  reference  point  y  called  the  “mother”  and  the  remaining  j 

spacecraft  denoted  as  “daughters”.  A  circular  formation  can  be  defined  even  in  the 

absence  of  a  mother  spacecraft.  Thus,  the  mother  spacecraft  may  be  replaced  by  a 

3 


reference  point,  y,  which  serves  the  purpose  of  providing  a  non-inertial  reference  frame 
to  the  entire  fonnation.  Multiple  rings  of  circular  formations  can  similarly  be  defined 
with  multiple  distance  bounds  on  a  collection  of  spacecraft  with  respect  to  a  single 
reference  point. 

A  fonnation  is  defined  to  be  fully  periodic  if  the  entire  state  vector  is  periodic, 

x'(x0)  =  x'(x/)  (2.6) 

and  partially  periodic  if  only  some  of  the  components  of  the  state  vector  are  periodic. 
For  example,  a  formation  may  be  partially  periodic  because  it  is  periodic  in  position  but 
not  in  velocity,  or  vice  versa.  If  propellant  is  used  to  maintain  a  formation,  then  by  this 
definition,  the  formation  is  partially  periodic  if  the  aperiodic  mass  is  included  as  a  state 
variable.  It  is  apparent  that  a  circular  formation  is  a  periodic  fonnation  but  the  reverse  is 
not  necessarily  true.  Of  special  note  is  that  this  definition  for  periodic  motion  can  be 
either  inertial  or  relative.  Hence,  periodicity  in  the  relative  frame  does  not  necessarily 
imply  periodicity  in  the  inertial  frame.  That  is,  the  swann  may  drift  in  inertial  space,  but 
it  will  stay  cohesive  as  a  formation. 

For  a  fully  periodic  fonnation  using  Cartesian  position  and  velocity  vectors  as  the 
state,  the  periodicity  constraint  may  be  written  as, 

r(x0)  =  r(x/)  (2.7) 

f(x0)  =  f(x/)  (2.8) 

These  conditions  allows  us  to  further  define  two  classes  of  partially  periodic  problems: 

(1)  when  only  equation  (2.7)  is  imposed  while  the  boundary  conditions  on  r  are  free  and 

(2)  when  only  equation  (2.8)  is  imposed  while  the  boundary  conditions  on  r  are  free. 
This  thesis  will  limit  its  attention  to  fully  periodic  solutions  and  the  two  classes  of 
partially  periodic  fonnations  described  above. 
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Finally,  it  must  be  noted  that  a  formation  need  not  be  periodic  at  all!  Hence,  a 
relaxed  formation  is  defined  as  the  case  when  the  state  vector  returns  to  within  a  defined 
space  around  the  initial  state  after  some  time  x/ , 


F  <  X1  -  x'  <  F 

b,  -  Aq  A j-  _  bu 


Thus,  the  familiar  notion  of  an  epsilon  ball, 


x‘  -x‘ 

0  / 


<£ 


(2.9) 


(2.10) 


is  included  in  this  definition. 


The  configuration  is  considered  to  be  optimal  if,  in  addition  to  satisfying  the 
configuration  constraint,  a  scalar  performance  measure, 


1  c 

J[x(-),u(-),x0,x/;p]=  — — —  JV(x(x),u(x),x;p)dx 


(2.11) 


is  minimal.  The  reason  for  choosing  a  cost  functional  borrowed  from  optimal  periodic 
control  theory  is  that  orbital  motion  is,  by  nature,  periodic.  Further,  in  addition  to  finding 

minimal  fuel  configurations,  it  is  also  desirable  to  find  the  optimal  period,  [x;  -x0J  that 

renders  a  minimal  cost  per  period. 

To  complete  the  optimal  periodic  control  formalism,14  the  configuration 
constraints  described  in  equation  (2.1)  are  broken  down  into  event  constraints  and  path 
constraints.  Event  constraints  are  end  point  boundaries  defined  by 

e/  - e (x1  (X0 ) , X1'  (x7 ),x0,x/)  < e„  (2. 12) 


Path  constraints  are  boundaries  placed  on  the  trajectory  of  the  model, 

g,  <g(x'  (x),u'  (x),x)<g„  (2.13) 


Additionally,  each  of  the  state  and  control  variables  may  have  a  constraint  placed  on 
them  by 
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(2.14) 


x*  <  x(  (x)  <  x1 

u  <u  (x)<u  (2.15) 

All  of  the  constraints  shown  above  can  be  used  as  equality  constraints  by  simply 
setting  the  upper  and  lower  bounds  equal.  They  are  written  as  inequalities  for  the 
purpose  of  generality.  Any  formation  configuration  may  now  be  defined  in  this 

“standard”  fonn  as  finding  the  controls,  u'(x),  and  the  optimal  period,  [x/-x0],  that 

minimize  the  cost  of  equation  (2.11). 

1.  Reference  Frame 

In  order  to  describe  relative  position,  motion,  and  configurations,  the  Formation 
Reference  Frame  will  be  used  throughout  this  thesis.  Figure  II- 1  shows  this  reference 
frame,  which  is  defined  with  x  pointing  in  the  radial  direction,  y  pointing  perpendicular 

to  x  along  the  direction  of  motion,  and  z  completing  the  right-handed  coordinate 

20 

system.  This  reference  frame  is  often  called  RSW  or  Satellite  Coordinate  System. 


6 


The  assumption  of  the  coordinate  system  is  important  in  defining  any  problem. 
For  some  coordinate  systems,  the  fonnation  relative  motion  and  configuration  constraints 
are  intuitive.  For  example,  the  path  constraint  for  a  circular  fonnation  can  be  described 
in  the  formation  reference  frame  simply  using  the  relative  position 

/;2  +  r2  +  rz2  =  constant  (2.16) 

On  the  other  hand,  to  describe  this  constraint  in  the  inertial  frame  requires  a 
transfonnation  to  that  frame.  This  transfonnation  need  only  be  completed  once,  but  can 
be  very  labor  intensive. 

This  same  concern  must  be  applied  to  the  dynamic  equations  of  motion.  The 
complete,  non-linear  equations  can  be  readily  expressed,  with  no  assumptions,  in  the 
inertial  frame  by 

<2-17> 

R 

Conversely,  in  the  formation  reference  frame  a  transformation  is  again  required  and 
usually  involves  introducing  assumptions  and  linearizations. 

Equal  to  the  coordinate  system  in  importance  is  the  choice  of  the  variables  used  to 
describe  the  satellite  state.  One  set  of  variables  is  the  Cartesian  position  and  velocity 
vectors.  Other  sets  available  include  many  different  orbital  element  sets,  which  are 
especially  useful  if  the  coordinate  system  is  Earth  centered  and  inertial.  However,  if  the 
solution  calls  for  control  thrusting,  it  will  be  a  complex  translation  into  the  orbital 
elements.  For  these  and  other  reasons,  all  models  presented  in  this  thesis  utilize  Cartesian 
position  and  velocity  vectors  in  the  Formation  Reference  Frame  as  the  basic  spacecraft 
state.  Depending  on  the  model,  there  may  be  additional  variables  in  the  state,  but  there 
will  be  a  minimum  of  these  six. 

B.  SOLVING  OPTIMAL  CONTROL  PROBLEMS 

Until  recently,  solving  general  nonlinear  optimal  control  problems  was  an  arduous 
or  impossible  task.  The  theoretical  framework  for  solving  such  problems  is  the  Minimum 
Principle.  Numerical  methods  based  on  the  Minimum  Principle  are  known  as  indirect 
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methods.  While  solutions  obtained  from  indirect  methods  are  accurate  in  the  sense  that 
they  satisfy  the  necessary  conditions  of  optimality,  they  are  fundamentally  burdened  by 
numerical  sensitivities  as  noted  by  Kalman  over  four  decades  ago."  The  so-called 
indirect  multiple  shooting  methods  and  indirect  collocation  methods  overcome  this 
computational  instability  problem  but  at  the  expense  of  convergence:  good  guesses  on  the 
costate  time-history  are  necessary  to  successfully  solve  the  problem.  Over  the  last 
decade,  the  so-called  direct  methods  have  come  to  the  forefront.  These  methods  bypass 
the  Minimum  Principle  and  “directly”  solve  the  problem  in  various  ways.  Betts  provides 
an  excellent  review  of  this  approach.  ~  Early  direct  methods  were  plagued  by 
inaccuracies,  particularly  in  the  determination  of  the  controls.  More  recently,  major 
breakthroughs  in  higher-order  methods  and  large  sparse  numerical  methods  have  quickly 
narrowed  this  gap.  ~  One  particular  approach  is  to  use  a  solution  obtained  from  a  direct 
method  as  a  guess  for  an  indirect  method.  Another  approach,  favored  in  this  thesis,  is 
called  the  Legendre  pseudospectral  method.  ~  This  method  is  used  to  solve  the  formation 
design  and  control  problems  posed  and  is  implemented  in  the  reusable  software  package, 
DIDO.  Unless  otherwise  specified,  all  results  reported  in  this  thesis  are  obtained  using 
this  software. 

1.  Necessary  Conditions  for  Optimality* 

The  first  step  in  solving  an  optimal  control  problem  is  to  construct  a  scalar 
function  called  the  Hamiltonian,  H , 

,  .  Efix,^)  .  ,T  /  , 

Ef(x,u,t,k)  = — - - -  +  k(7)  f(x,u,t)  (2.18) 

Tf  —  T0 

where  f(x,u,t)  are  the  dynamic  constraints  on  the  system,  and  k  ( l )  are  the  Lagrange 

multipliers  called  costates.  According  to  the  Minimum  Principle,  at  each  instant  of  time, 
the  optimal  control  is  obtained  by  solving  the  following  problem. 

Minimize  H  with  respect  to  u ,  with  u  e  U 


*  Most  of  the  information  in  this  section  comes  from  class  notes  and  discussion  from  AA  4850  and  is 
reproduced  here  for  completeness. 
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where  XJ  is  the  set  of  all  allowable  controls.  To  solve  this  problem,  the  Lagrangian  of 
the  Hamiltonian  must  be  constructed: 

L(x,u,fik,(|>)  =  H(x,u,  t,  X)  +  <>  (t)Tg(x,u,t)  (2.19) 

where  g(x,u,t)  are  the  path  constraints  and  (j)(/)  are  the  associated  Lagrange 

multipliers.  The  path  constraints  include  all  trajectory  path  constraints  as  well  as  any 
state  and  control  bounds.  Applying  the  Karush-Kuhn-Tucker  (KKT)  theorem  to  the 
Lagrangian  results  in  a  set  of  necessary  conditions  and  provides  a  means  to  demonstrate 
the  optimality  of  a  solution. 


^  =  0 

du 

(2.20) 

(j)(t)T  g(x,u,t)  =  0 

(2.21) 

with 

<0 

g,  =g(x,U,l) 

^  =  < 

“°  if 
=  0 

g  (XjU,  x)  =  gu 

term  bv  term 

g,  <g(x,u,x)<g„ 

(2.22) 

any 

g,  =  gu 

The  third  case  above  describes  a 

special  condition  when  the  constraints  in  g 

are  interior 

or  non-binding, 

g/<g(x,u,t)<g„ 

(2.23) 

For  these  cases,  the  multipliers,  (|) 

=  0  and  equation  (2.20)  simplifies  to 

dL  _  _  o 

du  du 

(2.24) 

It  is  desirable  to  have  interior  constraints  because  the  problem  behaves  as  if  the 
constraints  do  not  exist.  For  this  reason,  the  constraints  placed  on  a  problem  may  have  an 
actual  value  in  practice,  but  if  sufficiently  large  as  to  be  non-binding,  they  can  be 
described  as  unconstrained. 
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2.  Scaling 

Scaling  is  critical  to  all  optimal  control  problems.  The  goal  is  to  establish  a 
scheme  that  scales  ah  parameters  and  state  variables  so  that  they  are  close  to  one  another 
in  value.  For  example,  it  is  better  to  have  positions  from  0-10  with  velocities  from  1-10 
than  positions  from  7000-17000  with  velocities  from  0-10.  The  method  for  scaling  the 
problem  starts  with  setting  a  desired  unit  and  detennining  the  remaining  units  that 
comply  with  this  standard.  For  the  problems  that  follow,  the  desired  standard  unit  was 
time.  The  time  unit  was  chosen  to  go  from  0  to  1  for  a  single  orbital  period.  The 
nonnalization  constant,  TU ,  is  defined  to  be  equal  to  the  period  of  the  reference  orbit  in 
seconds.  This  can  be  done  several  ways,  but  the  most  simple  way  is  to  select  a  value  for 
the  semi-major  axis,  a  ,  and  use  the  following  equation 


TU  =  27t 


(2.25) 


where  |l  is  the  gravitational  parameter  of  the  earth  and  is  set  to 


|l  =  3.986004418xl014 

s 

The  normalized  time  now  becomes 


(2.26) 


x 


nondim 


x 

TU 


(2.27) 


and  therefore,  the  orbital  period  becomes  1.0  TU.  This  definition  for  the  time  unit,  by 
nature,  defines  the  mean  motion  since 


n 


nondim 


2n 

Orbital  Period 


(2.28) 


The  next  variable  to  be  selected  is  the  mass  unit,  MU ,  which  is  simply  chosen  to 
be  equal  to  the  initial  mass  of  the  spacecraft  or 


mn 


--^  =  1.0 


MU 


(2.29) 
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The  last  variable  that  is  chosen  by  the  user  is  the  standard  distance  unit,  RU  .  The 
choice  for  this  unit  is  based  on  the  desired  formation  spacing.  It  can  be  varied  to  provide 
larger  or  smaller  formations  without  affecting  the  configuration.  Given  TU  and  RU , 
the  remaining  variables  can  be  nonnalized  as  follows 


nondim 


r 

RU 


(2.30) 


1  nondim 


(2.31) 


a 


P  nondim 


(2.32) 


Unless  dimensional  units  are  specified,  all  values  are  assumed  to  be  scaled  and 
nondimensional  and  the  nondim  subscript  will  be  dropped. 

3.  Unique  Issues  for  Periodic  Problems 

As  with  any  optimization  problem,  most  errors  are  not  in  the  solution  but  in 
asking  the  right  question.  One  special  issue  with  periodic  optimization  problems  is 
formulating  an  initial  guess.  Every  optimizer  requires  a  first  guess,  whether  it  is  provided 
by  the  user  or  calculated  by  the  solver.  The  quality  of  the  guess  is  often  an  issue,  as  some 
solvers  require  relatively  accurate  guess.  DIDO  does  not  require  a  good  quality  guess  or 
even  a  feasible  guess.  However,  a  simple  linear  interpolation  between  the  endpoints  is 
not  an  option  for  periodic  problems.  The  explanation  lies  in  the  periodicity  condition 
itself.  For  orbits,  even  a  straight  line  cannot  be  used  as  a  guess  since  the  motion  is 
elliptical  in  nature.  Some  sort  of  ellipse  must  be  used  as  a  guess  for  orbital  motion. 

Determining  the  cost  function  is  another  difficult  part  of  configuration  design.  It 
affects  not  only  the  speed  of  the  solution,  but  the  solution  itself.  For  example,  beginning 
with  section  III.D,  thrust  and  mass  are  used  in  the  non-linear  equations  of  motion.  The 
cost  function  can  be  written  simply  as 


J 


mf  —m0 
T/-T» 


(2.33) 
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On  the  other  hand,  it  can  also  be  written  as 


X  ,  -  X 

J  °  t0 


(2.34) 


which  also  results  in  a  minimum  mass  expended.  It  turns  out  that  for  zero  or  very  low 
thrust  problems  the  change  in  mass  is  very  small  and  may  be  below  the  numeric  tolerance 
of  the  solver.  This  often  results  in  a  sub-optimal  solution.  When  thrust  is  used  in  the 
cost,  since  it  is  a  control  variable,  the  solver  is  better  able  to  find  the  optimal  solution. 

Another  aspect  that  proved  interesting  was  the  (x/-x0)  term,  used  to  represent 

the  optimal  period.  For  Natural  or  Zero-Thrust  solutions,  the  optimal  period  was  exactly 
equal  to  the  orbital  period.  However,  once  thrust  is  required  to  maintain  the  formation 
the  optimal  period  may  or  may  not  equal  the  orbital  period.  One  way  to  get  around  this 
disparity  is  to  force  the  solution  period  to  be  equal  to  the  orbital  period,  that  is  t  f  =  1.0  . 


C.  VALIDATING  SOLUTIONS 

The  purpose  of  this  section  is  to  describe  the  various  methods  used  to  validate  the 
numerical  solutions.  Numerical  methods  may  seem  to  discard  the  physics  of  the  problem 
and  rely  solely  on  the  mathematics,  but  this  is  far  from  the  truth.  Any  solution  found 
numerically  must  be  ‘assessed’  to  see  if  it  is  in  fact  a  feasible  and  legitimate  solution.  It 
is  precisely  here  that  the  “Physics”  intuition  and  knowledge  are  implemented  to  assist  in 
verification.  To  that  end,  there  are  many  different  ways  to  verify  that  a  solution  is 
feasible.  The  optimality  of  a  solution  can  be  demonstrated  using  the  optimality 
conditions  described  in  section  B.l.  Each  of  the  three  primary  means  of  validating 
solutions  are  described  below. 

1,  Numeric  Propagators 

The  primary  method  for  validating  feasibility  of  the  solutions  is  to  propagate  them 
using  a  numerical  propagator.  Since  the  equations  of  motion  are  Ordinary  Differential 
Equations  (ODE),  a  tool  is  needed  to  solve  them  for  each  desired  time  step.  The 
environment  for  DIDO  is  MATLAB8,  so  a  resident  ODE  solver  was  used.  There  are 
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several  to  choose  from,  but  ODE45  is  the  most  common  and  usually  the  most  capable  for 
non-stiff  equations. 

The  first  step  in  validating  is  to  build  a  propagator  from  the  ODE  solver.  This 
requires  defining  the  Equations  of  Motion  (EOM)  that  govern  the  behavior  for  the  solver 
and  establishing  the  initial  conditions.  The  initial  conditions  were  taken  directly  from  the 
DIDO  solution.  It  is  important  for  the  EOM  to  be  the  same  as  the  ones  used  to  find  the 
solution,  to  ensure  an  accurate  validation  under  the  same  assumptions.  It  does  no  good  to 
make  assumptions  in  defining  the  problem,  only  to  forget  about  them  in  validating  the 
results. 

2.  Commercial  Software 

Another  method  employed  to  validate  solutions  is  the  use  of  commercial  software. 
Both  Satellite  Tool  Kit®  (STK),  from  Analytical  Graphics,  Inc.  and  FreeFlyer®,  from  A. I. 
Solutions  were  used  as  an  independent  source  for  propagating  the  formations.  Usually 
these  programs  are  used  for  visualization  and  presentation,  but  both  software  suites 
include  a  very  robust  non-linear  propagator.  Any  number  of  reports  can  be  selected  to 
identify  and  track  the  desired  parameters.  Each  of  the  programs  also  includes  a  viewing 
capability  that  allows  the  user  to  be  placed  on  the  reference  point  observing  the  formation 
from  that  perspective.  This  directly  corresponds  to  the  Formation  Reference  Frame 
shown  in  Figure  II- 1  and  is  very  useful  in  visualizing  relative  fonnations. 

The  mechanism  for  validating  solutions  using  these  programs  is  to  import  the 
initial  conditions  for  the  individual  fonnation  satellites  into  the  programs  and  propagate 
them  forward.  This  can  be  done  directly  from  MATLAB  for  use  in  FreeFlyer,  or  with  an 
Ephemeris  file  for  use  in  STK.  If  the  solution  requires  active  control,  then  the  control 
history  must  also  be  imported  into  the  programs.  Other  than  initial  conditions  and/or 
controls,  there  is  no  other  information  provided  to  the  software.  This  lack  of  information 
is  exactly  what  validates  the  solutions.  If  the  formation  behaves  as  predicted  in  the 
solution  when  propagated  by  the  software,  then  the  formation  is  feasible. 
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3.  Previously  Discovered  Solutions 

The  third  method  is  used  to  validate  the  method  more  than  it  is  to  validate  a 
particular  solution.  If  the  method  is  employed  and  is  able  to  identify  previously 
discovered  solutions,  it  is  valid.  The  ultimate  validation  would  be  to  demonstrate  that 
this  method  is  able  to  find  all  previously  discovered  solutions.  Rather  than  devote 
precious  time  to  an  exhaustive  catalog  of  all  prior  solutions  to  all  satellite  formations, 
several  representative  formations  were  chosen  as  demonstrations  of  the  ability. 
Specifically,  the  well-known  circular  and  projected  circular  solutions  to  the  Hill- 
Clohessy-Wiltshire  equations  were  reproduced.  For  an  example  using  elliptical  reference 
orbits,  Reference  17  outlines  a  solution  with  e  =  0.7  for  the  reference  satellite. 
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III.  CIRCULAR  REFERENCE  ORBIT 


A.  PROBLEM  DEFINITION 

The  framework  described  in  the  previous  chapter  can  be  applied  to  designing  and 
controlling  spacecraft  fonnations  subject  to  any  arbitrary  forces.  This  chapter  begins  by 
applying  this  framework  to  spacecraft  subject  to  an  inverse-square  gravity  field  only.  If 
the  equations  of  motion  are  linearized  in  the  formation  reference  frame  about  a  circular 
reference  orbit,  the  C-W  equations  are  obtained  and  closed-form  solutions  are  easily 
found.  From  these  equations,  it  is  apparent  that  zero-propellant  fonnations  exist  making 
this  problem  an  excellent  starting  point.  These  equations  also  served  as  a  good  model 
with  which  to  validate  the  process  for  this  class  of  periodic  problems. 

One  zero-propellant  solution  to  the  C-W  equations  is  a  circular  formation.  In 
order  to  solve  this  problem  the  following  configuration  was  used.  The  state  consists  of 
the  Cartesian  relative  position  and  velocity  components  shown  below 

x  =  [rT  i*T]T  =[r.  r,  rx  rx  ry  r,J  (3.1) 

The  state  constraints  were 


-2 r < r  r  r  <  2  r 


-r  <  r  r  r  <  r 

max  x  ’  y  ’  z  max 


(3.2) 


where  rmax  was  chosen  arbitrarily.  The  position  constraints,  while  specified,  were  never 
active  due  to  a  tighter  path  constraint.  The  velocity  constraints  were  not  active  either  due 
to  the  choice  of  rmax .  The  controls,  representing  the  relative  accelerations  in  each  axis, 
were  constrained  by 


— u  <u  ,u  ,u  <u 

max  —  x  ’  v  ’  z  —  n 


(3.3) 


with  umax  specified  by  the  user.  Since  the  solution  is  a  zero  control  solution,  these 
constraints  were  also  inactive. 
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The  event  constraints  were  specified  for  a  fully  periodic  solution  using  equations 
(2.7)  and  (2.8)  although  they  were  written  in  the  form 

x(x0)-x(x/)  =  0  (3.4) 

or  more  specifically 

f(x/) 

rv(X/) 

fi(x/) 
fi(x/) 
fi(x/) 
f(x/) 

The  path  constraint,  g  ,  for  a  circular  fonnation  is  defined  by 

(r-8)2  <(rc2+ry2  +  rz2)<(r  +  8)2  (3.6) 

This  general  form  of  the  path  constraint  includes  an  allowance  for  a  tolerance,  which  for 
circular  formations  is  set  to  zero  or  for  near-circular  formations  is  nonzero.  For  this  case, 
the  path  constraint  for  a  circular  formation  and  the  event  constraint  for  a  fully  periodic 
formation  are  redundant.  The  path  constraint  will  force  a  fully  periodic  solution,  even  in 
the  absence  of  any  event  constraints.  This  was  demonstrated  and  validated  during  several 
of  the  solution  sets. 

The  cost  function  to  be  minimized  was 

1  v 

J  = - f(w"+w~+zcWx  (3.7) 

X  ,  —  X  J  V  '  ' 

J  o  x0 

which  is  identical  to  a  cost  that  was  based  on  the  absolute  value  of  each  control 
acceleration,  if  the  solution  is  a  zero-control  solution. 


"fi(Xo)" 
ry  (X«) 

f(Xo) 

fi(xo) 

f(xo) 

_fi(xo)_ 
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B.  EQUATIONS  OF  MOTION 

Since  the  equations  of  motion  for  this  model  are  well  known1 5,24,25,26  their 
derivation  will  not  be  described  here.  Instead,  only  the  equations  and  their  assumptions 
will  be  presented. 


~  0 

2  n 

O' 

"3  nr 

0 

0  " 

u. 

r  = 

-2  n 

0 

0 

r  + 

0 

0 

0 

r  + 

Uy 

0 

0 

0 

0 

0 

2 

—n 

(3.8) 


where  r  is  the  vector  representing  the  inter-satellite  distance  expressed  in  the  formation 
reference  frame  (see  Figure  II- 1), 


r  =  [x  y  z  f 

or  in  state  form: 


"  0 

0 

0 

1 

0 

0  " 

"0" 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

1 

x  + 

0 

3  n2 

0 

0 

0 

2  n 

0 

0 

0 

0 

-2  n 

0 

0 

Uy 

_  0 

0 

-«2 

0 

0 

0 

(3.9) 


(3.10) 


Two  of  the  primary  assumptions  of  this  model  are  a  spherical  earth  with  no  other 
perturbing  forces  present  and  a  circular  reference  orbit.  Both  of  these  assumptions  lead 
to  unstable  fonnations  when  applied  to  actual  orbits  since  their  effects  are  significant. 
The  latter  assumption  is  the  primary  focus  of  this  thesis  and  will  be  addressed  in  a  later 
chapter.  The  third  assumption  is 

r«R  (3.11) 


which  remains  valid  for  most  formations  even  when  applied  to  actual  orbits. 

C.  C-W  SOLUTIONS 

For  all  solutions  shown  in  this  thesis,  the  filled  circles  show  the  node  points 
corresponding  to  the  solution.  They  vary  in  number  based  on  an  arbitrary  user-defined 

specification.  The  nature  of  the  solver  is  such  that  the  spacing  between  node  points  is  not 
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constant.  Instead,  the  points  are  closer  together  at  the  endpoints  than  they  are  in  the 
middle.  The  solid  lines  represent  the  motion  of  the  propagated  initial  conditions.  The 
solutions  for  the  C-W  equations  of  motion  are  not  shown  for  the  sake  of  brevity  since 
they  are  identical,  within  numeric  tolerances,  to  the  solutions  for  the  new  equations 
shown  later. 

1.  C-W  Circular  Solution 

For  the  C-W  equations,  imposing  the  condition  of  full  periodicity  and  restricting 
the  inter-satellite  range  to  a  constant  distance  from  the  center  X  (at  point  [0,0,0]) 
produces  a  circular  satellite  formation. 

The  final  time  for  the  optimal  period  was  not  fixed.  The  solution  was  not  only 
able  to  detennine  the  optimal  trajectory,  but  also  the  optimal  period  in  which  to  complete 
its  trajectory.  As  expected,  the  solution  resulted  in  xf  =  1.0 ,  meaning  the  optimal  period 

was  exactly  equal  to  one  orbital  period.  Interestingly,  when  asked  to  find  a  solution  over 
2  orbits,  the  same  1 -orbit  solution  was  found,  only  it  was  now  shown  over  the  2  orbits. 

2.  C-W  Projected  Circular  Solution 

Another  well-known  solution  to  the  C-W  equations  is  the  projected  circular 
formation.  This  solution  maintains  the  appearance  of  a  circular  formation  as  seen  from 
the  surface  of  the  earth,  but  is  in  fact  elliptical.  To  accomplish  this  fonnation  a  different 
path  constraint,  g  ,  was  imposed. 

(r  — 8)”  <  +r.2)  <  (r  +  S)2  (3.12) 

It  is  no  longer  the  three  dimensional  range  that  was  constrained  but  the  range  from 
the  reference  point  to  the  satellite  in  a  certain  plane,  namely  the  cross-track  versus  along- 
track  plane.  As  with  the  circular  fonnation,  the  optimal  period  was  not  fixed,  and  again  it 
was  exactly  equal  to  the  orbital  period. 

D.  NON-LINEAR  TWIST 

To  further  amplify  the  notion  that  linear  models  are  not  necessary  for  this 
approach,  Thrust  was  chosen  as  a  control  variable  instead  of  acceleration.  Not  only  is 
this  more  realistic  but  it  also  makes  the  “linear”  equations  “nonlinear”  due  to  a  non-zero 
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mass  flow  rate.  This  requires  the  addition  of  Mass  as  one  of  the  state  variables,  in 
addition  to  the  relative  position  and  velocity  components. 


x  = 


(3.13) 


The  path  constraint  remains  the  same  as  equations  (3.6)  or  (3.12).  The  state 
constraints  now  include  the  following  constraint  on  mass: 

0.10  <m(x)<  1.0  (3.14) 


which  assumes  a  minimum  mass  of  10  percent  of  the  original  mass  and  a  maximum  mass 
equal  to  the  original  mass.  The  event  constraints  also  contain  the  new  event 

m(x0)  =  1.0  (3.15) 

to  force  a  starting  mass  at  the  beginning  of  the  solution. 

The  controls  represent  thrust  in  a  given  axis  and  are  constrained  by 

-T^T„T,,Tz<T„  (3. 16) 


with  Tmax  usually  chosen  in  such  a  way  as  to  make  the  constraint  inactive.  The  desired 
cost  is  the  total  amount  of  thrust  expended  and  can  be  described  as 

1  v 

J  = - \\Tx\  +  \T\  +  \Tz\dx  (3.17) 

x/-x»  i  1  1 

but  was  implemented  as 

1  T/ 

J  = - J  (T;  +T;  +T;)dT  (3.18) 

T/  ~  X0  To 

This  cost  function  could  have  been  implemented  as  seen  in  equation  (3.17)  but,  for 
numerical  performance  reasons,  equation  (3.18)  was  chosen.  If  the  cost  in  equation 
(3.17)  is  minimized,  the  fonnation  solution  will  be  identical  to  the  formation  found  by 
minimizing  the  cost  in  equation  (3.18),  when  7)  =  0  .  For  zero-thrust  solutions,  the  main 
difference  is  the  speed  of  the  process. 
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1.  New  Equations  of  Motion 

The  equations  of  motion  for  this  case  are  similar  to  Equation  (3.8)  except  the 
control  accelerations,  u ,  are  replaced  with  thrust. 


~  0 

2  n 

O' 

~3n2 

0 

0  " 

[T~\ 

1 

X 

r  = 

-2  n 

0 

0 

+ 

•u 

0 

0 

0 

r  H — 

m 

Ty 

0 

0 

0 

0 

0 

-n1 

Jz_ 

(3.19) 


Another  equation  must  be  included  to  govern  the  new  state  variable  of  mass.  This  mass 
flow  rate  is 


-T 


(3.20) 


Note  that  T  is  not  the  magnitude  of  the  thrust  vector,  but  is  defined  as 
T  =  T  |  +  |r  |  +  It  |  and  ve  is  the  characteristic  exhaust  velocity  of  the  thruster. 


2.  Circular  Formation  Solution  with  New  EOM 

The  circular  formation  solutions  to  the  new  EOM  are  obtained  using  the  path 
constraint  shown  in  equation  (3.6).  Table  III- 1  shows  the  constraints  used  for  this 
solution.  Figure  III- 1  shows  the  three-dimensional  view,  while  Figure  III-2  to  Figure 
III-4  show  the  orthogonal  projections.  Figure  III-3  also  shows  the  planar  motion  in  the 
x  —  z  plane  at  an  angle  of  60°  to  greater  than  10'9  accuracy.  Figure  III-5  shows  the  thrust 
profile,  which  is  equal  to  zero-thrust  within  numerical  tolerances. 
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Cross-Track 


Table  III- 1  Constraints  for  Circular  Fonnation 


Constraint 

Normalized 

Lower  and  Upper  Bounds 

States:  rx,ry,rz 

Unconstrained 

States:  rx,ry,rz 

Unconstrained 

States:  m 

[0.1  :  1.0] 

Controls:  T 

Unconstrained 

Time:  x 

Unconstrained 

Events:  r(x0)-r(x/-) 

[0.0  :  0.0] 

Events:  r(x0)-r(x/-) 

[0.0  :  0.0] 

Path:  g  =  r  +  r2  +  r 2 

[1.0  :  1.0] 

Number  of  Nodes 

120 

Reference  Orbit:  e 

0.0 

Figure  III- 1  Circular  Formation  Using  New  Equations  of  Motion 
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Cross-Track 


-0.8  -0.6  -0.4  -0.2  0  0.2  0.4  0.6  0.8 

Along-Track 


Figure  III-2  Radial  vs.  Along-Track  Motion  for  a  Circular  Formation 


Radial 


Figure  III-3  Cross-Track  vs.  Radial  Motion  for  a  Circular  Fonnation 


Figure  III-5  Thrust  Profile  for  Circular  Formation  with  New  EOM 
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3.  Projected  Circular  Formation  Solution  with  New  EOM 

The  projected  circular  formation  solutions  utilize  the  new  EOM  with  the  path 
constraint  in  equation  (3.12).  Table  III-2  shows  the  constraints  while  Figure  III-6  shows 
the  three-dimensional  view.  Figure  III-7  to  Figure  III-9  show  the  orthogonal  projections 
and  Figure  III- 10  shows  the  thrust  profile. 


Table  III-2  Constraints  for  Projected  Circular  Formation 


Constraint 

Normalized 

Lower  and  Upper  Bounds 

States:  r  ,r  ,r 

Unconstrained 

States: 

Unconstrained 

States:  m 

[0.1  :  1.0] 

Controls:  T 

Unconstrained 

Time:  x 

Unconstrained 

Events:  r(x0)-r(x/.) 

[0.0  :  0.0] 

Events:  r(x0)-r(x/.) 

[0.0  :  0.0] 

Path:  g  =  r2  +  r2 

[1.0  :  1.0] 

Number  of  Nodes 

100 

Reference  Orbit:  e 

0.0 
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Radial 


Figure  III-6  Projected  Circular  Formation  Using  New  Equations  of  Motion 


Figure  III-7  Radial  vs.  Along-Track  Motion  for  a  Projected  Circular  Formation 


Figure  III- 1 0  Thrust  Profile  for  Projected  Circular  Formation  with  New  EOM 


E.  VALIDATION  OF  SOLUTIONS 

Since  the  answer  to  these  problems  was  known,  the  validation  was  simply  to 
compare  the  numerical  solution  to  the  analytical  solution.  This  was  accomplished  using 
Figure  III- 1  through  Figure  III- 10  as  well  as  unique  characteristics  such  as  the  angle  of 
the  plane  of  relative  motion  seen  in  Figure  III-3.  This  angle  was  calculated  for  each 
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solution  and  verified  against  the  known  values." 

Validation  was  also  accomplished  by  numerically  propagating  the  initial 
conditions  for  a  set  number  of  orbits,  arbitrarily  chosen  as  50  orbits  and  is  shown  in 
Figure  III- 1 1  and  Figure  III- 12.  Table  III-3  shows  the  results  from  propagating  the  initial 
conditions  of  the  analytical  solution.  Table  III-4  and  Table  III-5  show  the  results  for 
propagating  the  initial  conditions  from  the  DIDO  solution  for  the  circular  formation  and 
projected  circular  formation,  respectively.  The  errors  seen  in  Table  III-3  fonn  the 
baseline  for  interpreting  the  errors  in  every  other  solution  propagation.  The  source  of 
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these  errors  is  the  numeric  propagator,  since  the  equations  of  motion  used  to  propagate 
are  the  same  as  the  analytical  model. 

As  mentioned  previously,  the  optimal  period  was  not  fixed  for  these  problems. 
Each  of  the  solutions  in  this  chapter  yielded  an  optimal  period  equal  to  the  orbital  period, 
which  was  further  verification  of  agreement  with  the  known  solutions. 
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Cross-Track 


Figure  III- 12  Projected  Circular  Solution  Over  50  Orbits 


Table  III-3  Propagation  Results  for  Analytical  C-W  Circular  Solution 


State  Variable 

%  Difference  Between 
Initial  and  Final  Values 

C 

9.84  xlO'4 

ry 

1.33  xlO'5 

rz 

9.84  xlO'4 

K 

1.33  xlO'5 

c 

9.84  xlO'4 

K 

1.33  xlO'5 
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Table  III-4  Propagation  Results  for  Circular  Formation 


State  Variable 

%  Difference  Between 
Initial  and  Final  Values 

C 

9.32  xlO'4 

ry 

3.15  xlO'4 

rz 

9.32  xlO'4 

K 

3.15  xlO'4 

c 

9.32  xlO'4 

3.15  xlO'4 

Table  III-5  Propagation  Results  for  Projected  Circular  Formation 


State  Variable 

%  Difference  Between 
Initial  and  Final  Values 

C 

9.30  xlO'4 

c 

3.13  xlO'4 

C 

9.30  xlO'4 

K 

3.16  xlO'4 

c 

9.30  xlO'4 

c 

3.16  xlO'4 
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IV.  ELLIPTIC  REFERENCE  ORBIT 


A.  PROBLEM  DEFINITION 

In  many  applications,  it  is  desirable  to  design  fonnations  with  non-circular 
reference  orbits.  The  equations  of  motion  from  Chapter  III  are  inapplicable  for  elliptical 
reference  orbits.17'27  This  chapter  addresses  these  fonnations;  in  particular,  the  design  of 
new  formations  using  the  method  described  in  Chapter  II. 

Using  the  same  coordinate  system  as  the  previous  chapter,  an  additional  state  was 
added  to  simplify  the  equations  of  motion.  Adding  v  ,  or  True  Anomaly,  the  state  vector 
becomes 


The  controls,  representing  both  a  positive  and  negative  thrust  direction  for  each 
axis,  are 

»=[>„  T'_  Tn  Tr_  TZ_J  (4.2) 


and  were  constrained  by, 


0  <TX+,TX. 


T  T  T 
’U+’U-’  z 


z+’ 


T  <  T 


(4.3) 


The  rationale  for  choosing  the  controls  above  was  mainly  for  mathematical  and 
numerical  purposes.  In  calculating  both  the  mass  flow  and  cost  of  the  solution,  it  is 
necessary  to  calculate  a  total  thrust  value.  Normally  this  is  done  by  taking  the  sum  of  the 
absolute  value  for  the  individual  thrust  variable  which  presents  mathematical  challenges 
due  to  non-differentiability.  That  is,  the  absolute  value  of  any  variable  is  defined  as 


x\ 


A 


X 

-x 


x>0 
x  <  0 


(4.4) 


It  is  evident  that  at  x  =  0  ,  the  absolute  value  function  becomes  nondifferentiable. 
This  can  be  avoided  by  using 


x\ 
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(4.5) 


but  even  this  introduces  numerical  issues  since  the  derivative  of  the  square  root  function 
is  undefined  at  x  =  0 ,  which  is  of  great  concern  since  Ti  =  0  is  not  only  possible,  but 
desired. 

One  workable  solution  is  to  use  a  different  thrust  variable  for  the  positive  and 
negative  directions.  This  choice  models  real  life  more  closely  since  actual  thrusters  are 
capable  of  firing  only  in  one  direction.  More  importantly,  it  removes  the  numerical 
difficulties  by  allowing  total  thrust  to  be  defined  as 

T  =  I  »'  =  Tx+  +  Tx_ _  +  Ty+  +  Ty_  +  L+  +  Tz_  (4.6) 

This  implementation  does  produce  an  additional  concern  since  Tf  =  0  is  no  longer 

internal  to  the  constraints  on  the  controls  therefore  making  the  constraints  active. 
Activation  of  the  Lagrange  multiplier  for  the  controls  does  not  change  the  problem,  but  it 
does  make  the  analysis  and  validation  a  bit  more  involved.  Note  that  T  does  not 
represent  the  magnitude  of  the  thrust  vector.  It  has  been  redefined  according  to 
equation  (4.6). 

For  this  formation,  the  path  constraint  is  defined  by, 

dj  <  ( r 2  ( X)  +  r 2  ( X)  +  r.2 ( X) )  <  d2u  Vx  (4.7) 

where,  di  and  du  are  minimum  and  maximum  distance  from  the  reference  point  to  any 
spacecraft  in  the  swarm.  The  cost  function  is 

1  Xf 

J  = - \Tdx  (4.8) 

v-v. 

It  is  obvious  that  if  the  optimal  cost  turns  out  to  be  zero,  then  the  solution  corresponds  to 
a  zero-propellant  fonnation;  otherwise,  the  optimal  (i.e.  minimum  fuel)  open-loop 
controls  to  achieve  the  desired  formation  are  obtained. 

B.  EQUATIONS  OF  MOTION 

Since  dynamics  from  Chapter  III  are  invalid  for  an  elliptical  reference  orbit,  a  new 
set  of  equations  must  be  used.  These  equations  of  motion  are  described  in  References  17 
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and  27  and  are  derived  here  for  the  purpose  of  completeness.  Following  Kane’s 
notation,  the  superscript  N  is  used  to  denote  the  Newtonian  or  Inertial  reference  frame. 
Starting  with  the  most  general  equation  describing  the  motion  of  the  unperturbed 
reference  orbit, 


K=~v 


R 

|R|3 


(4.9) 


The  motion  of  a  swarm  satellite  is  given  by 


R  +  r 


|3  P  T 


(4.10) 


where 


R  +  if  =  (|r|3  +2(R.r)|R|  +  |r|2  |r|)7 


(4.11) 


If  the  assumption  R  »  r  is  now  introduced, 


R  +  r  1 


|R  +  r|3  iRl3 


„  .  R*r 

R  +  r  —  3  — -R 

|R|' 


(4.12) 


The  relative  dynamics,  though  still  in  the  inertial  frame,  can  now  be  described  by 

f  \ 


—HN  -iiN  — 


R 


"-3^R 
|R| 


+  a>a  nt 


(4.13) 


J 


In  order  to  retrieve  the  relative  dynamics  in  the  formation  reference  frame,  the 
Transport  Theorem  for  moving  coordinate  systems  must  be  applied.  The  formation 
reference  frame  is  denoted  by  the  superscript  B  .  Generally  it  is 


r,v  =rB  +  (w(b5xr)  +  (A'(oBxA'e/xr)  +  2(A'c}j8xfii)  (4.14) 

Solving  for  rB  yields 

fB=fA,-(wc>'xr)-(A,(o'xA,to'xr)-2(A'oJ'xfB)  (4.15) 
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For  an  unperturbed  orbit  around  a  central  body,  using  the  perifocal  inertial  coordinate 
system. 


N  B 

(O 


[0  0  v]  and  N6)B  =  [0  0  v 


•  IT 


(4.16) 


The  above  equations  imply  that  there  cannot  be  a  satellite  at  the  reference  point,  since  it 
is  defined  as  an  unperturbed  reference  orbit.  The  definition  of  angular  momentum  for 

9Q 

any  unperturbed  orbit, 


h  -  R2 v  =  =  na2\[i 


-e 


solving  for  v , 


na 


v  =  — —  Vl  -  e2  with  a 

R2 


f?(l  +  ecos(v)) 

(i-!) 


finally  yielding 


«(l-ecos(v))~ 


(i -*r 

Differentiating  equation  (4.19)  with  respect  to  time  produces 

..  -2«esin(v)v 


v  ■ 


-(l  +  ecos(v)) 


One  additional  tenn,  aT ,  is  defined  as  follows 


(4.17) 


(4.18) 


(4.19) 


(4.20) 


—  — 


_L 

m 


r  t~ 

r  r. 

-T  1 

X 

_  l 

x+ 

x- 

T 

T  -T 

m 

y 

y- 

T 

T  —T  1 

L  z  J 

L  ^ 

Z-  J 

(4.21) 
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Substituting  equations  (4.13),  and  (4.16)  through  (4.21)  into  equation  (4.15) 
provides  the  relative  dynamical  equations  of  motion  in  the  relative  coordinate  frame. 
They  are 


with 


~  0 

2v 

O' 

V+2  k  v 

0  " 

[T~\ 

fB  = 

-2v 

0 

0 

rB  + 

-v  v2-k 
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k  =  n 


l  +  ecos(v) 

l-e2 


(4.23) 


and  an  additional  dynamic  equation  needed  for  the  mass  flow, 

m  =  —  (4.24) 


One  significant  fact  is  that  these  equations  are  non-linear  in  mass.  In  order  to 
limit  the  scope  of  this  chapter,  problems  were  confined  to  a  zero  perturbation  case  by 
setting  ap  =  0  .  It  should  be  noted,  however,  that  the  effects  of  any  modelable  perturbing 

force  may  be  included  in  a  as  desired  to  increase  the  fidelity  of  the  model. 


C.  NATURAL  SOLUTIONS 

The  first  set  of  solutions  are  Natural  solutions  for  varying  eccentricities  of  the 
reference  orbit.  Natural  is  meant  to  denote  that  they  require  no  thrust  to  maintain 
configuration  and  are  at  least  partially  periodic  satisfying  the  traditional  formation 
preconceptions.  They  are  of  course,  only  as  good  as  the  dynamical  model  used  to  find 
them. 


1.  Natural  Formation  for  e  =  0.5 

The  solution  was  found  in  nonnalized  units,  allowing  it  to  be  applied  for  any 
desired  formation  distance,  RU .  For  the  purpose  of  clarity  and  comparison,  the 
solutions  will  be  presented  in  nondimensional  units.  Table  IV- 1  itemizes  the  constraints 
in  place  for  this  fonnation. 
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Figure  IV- 1  shows  the  three-dimensional  view  of  the  fonnation  relative  motion. 
Figure  IV -2  shows  the  orthogonal  projection  of  the  relative  orbit  in  the  radial  versus 
along-track  plane.  Figure  IV-3  shows  the  projection  in  the  cross-track  versus  radial  plane 
while  Figure  IV-4  projects  the  orbit  onto  the  cross-track  versus  along-track  plane.  As 
with  the  previous  figures,  the  X  is  the  location  of  the  reference  point.  From  these  results, 
position  appears  to  be  naturally  constrained.  In  other  words,  the  relative  orbit  is  fully 
periodic  even  though  the  constraints  are  imposed  for  partial  periodicity  in  velocity.  This 
can  be  explained  by  the  fact  that  a  fully  periodic  solution  would  satisfy  the  constraints  for 
a  partially  periodic  solution.  Additionally,  since  velocity  was  constrained,  natural  orbital 
motion  will  tend  to  drive  position  toward  the  periodic  solution.  This  is  not  a  law,  since  it 
is  possible  to  find  a  formation  that  is  periodic  in  velocity  but  not  in  position.  On  the  other 
hand,  if  the  position  is  constrained  with  velocity  free,  the  formations  will  return  to  their 
original  position  but  the  velocity  may  be  so  large  as  to  make  this  return  impossible.  This 
is  of  course,  “undesirable”  for  fonnation  configurations. 


Table  IV-1  Constraints  for  Natural  Fonnation  with  e  =  0.5 


Constraint 

Normalized 

Lower  and  Upper  Bounds 

States:  r  ,r  ,r 

Unconstrained 

States: 

Unconstrained 

States:  m 

[0.1  :  1.0] 

Controls:  T 

[0  :  Unconstrained] 

Time:  x 

[0:5] 

Events:  r(x0)-r(x/) 

Unconstrained 

Events:  r(x0)-r(x/.) 

[0.0  :  0.0] 

Path:  g  =  r;  +  r]  +  rz2 

[2  :  22] 

Number  of  Nodes 

199 

Reference  Orbit:  e 

0.5 
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Figure  IV- 1  3 -Dimensional  Formation  Trajectory  for  e  =  0.5  Natural  Formation 
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Figure  IV -2  Radial  vs.  Along-Track  Motion  for  e  =  0.5  Natural  Formation 
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Figure  IV- 5  Thrust  Profile  for  e  =  0.5  Natural  Formation 

Of  course,  all  these  results  are  within  numerical  tolerances  arbitrarily  chosen  to  be 
10'  in  normalized  units.  It  turns  out  that  the  error  in  the  position  vector  is  within  10' 
units  after  one  orbit.  From  Figure  IV-5,  above,  it  is  clear  that  this  is  a  zero-propellant 
formation  configuration. 

2.  Natural  Formation  for  e  =  0.3 

Solution  Two  is  another  new  formation  configuration  using  the  constraints  shown 
in  Table  IV-2,  specifically  with  e  =  0.3 .  Similar  to  solution  one,  the  time  unit,  TU  ,  was 
calculated  for  a  reference  orbit  with  a  perigee  altitude  equal  to  1000  km.  One  of  the 
differences  between  this  fonnation  and  the  previous,  aside  from  reference  orbit 
eccentricity,  is  the  periodicity  constraints.  This  solution  was  constrained  to  be  periodic  in 
position  with  velocity  free.  The  previous  solution,  in  section  1,  was  periodic  in  velocity 
with  position  free.  Figure  IV-6  shows  the  three-dimensional  trajectory  for  this  orbit  in 
the  formation  reference  frame.  Figure  IV-7  to  Figure  IV-9  show  the  orthogonal 
projections  of  the  relative  orbit.  Again,  the  X  in  the  figures  is  the  location  of  the 
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reference  point.  Finally,  Figure  IV- 10  demonstrates  that  this  is  another  zero-propellant 
formation. 
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Cross-T  rack 


Table  IV-2  Constraints  for  Natural  Fonnation  with  e  =  0.3 


Constraint 

Normalized 

Lower  and  Upper  Bounds 

States:  r  ,r  ,r 

Unconstrained 

States:  rx,ry,rz 

Unconstrained 

States:  m 

[0.1  :  1.0] 

Controls:  T 

[0  :  Unconstrained] 

Time:  x 

[0:5] 

Events:  r(x0)-r(x/.) 

[0.0  :  0.0] 

Events:  f(x0)-r(x/) 

Unconstrained 

Path:  g 

Unconstrained 

Number  of  Nodes 

80 

Reference  Orbit:  e 

0.3 

Figure  IV-6  3-Dimensional  Formation  Trajectory  for  e  =  0.3  Natural  Formation 
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3.  Natural  Formation  for  e  =  0.7 

Recently,  Inhalan  et  al  presented  new  periodic  formations  by  an  analytic  method. 
Solution  Three,  shown  in  Figure  IV-11  and  Figure  IV-12,  reproduces  a  particular 
formation  given  in  Reference  17  for  e  =  0.7  using  this  method. 

By  finding  the  same  solution,  the  method  is  again  validated  against  independent 
results.  The  solution  was  found  by  constraining  the  position  states  to  a  three-dimensional 
box  slightly  larger  than  the  proposed  solution  and  using  the  fully  periodic  event 
constraints  (see  Table  IV-3). 


Table  IV-3  Constraints  for  Natural  Fonnation  with  e  =  0.7 


Constraint 

Normalized 

Lower  and  Upper  Bounds 

States:  r  ,r  ,r 

[±0.62,  ±1.4,  -5. 8:1.0] 

States: 

Unconstrained 

States:  m 

[0.1  :  1.0] 

Controls:  T 

[0  :  Unconstrained] 

Time:  x 

[0  :  10] 

Events:  r(x0)-r(x/) 

[0.0  :  0.0] 

Events:  f(x0)-r(x/) 

[0.0  :  0.0] 

Path:  g 

Unconstrained 

Number  of  Nodes 

99 

Reference  Orbit:  e 

0.7 
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D.  FORCED  SOLUTIONS 

This  section  provides  a  glimpse  into  the  power  of  this  method.  Specifically,  that 
it  can  find  not  only  natural  formations,  but  also  controlled  or  forced  fonnations.  These 
solutions  are  unique  in  the  sense  that  Keplerian  orbital  dynamics  cannot  produce  these 
formations  since  they  must  be  actively  controlled.  This  method  will  determine  the 
feasibility  of  the  formation  and  will  automatically  provide  the  open  loop  controls  required 
to  maintain  the  desired  configuration. 

1.  Forced  Circular  Formation 

The  fonnation  in  this  section  uses  a  reference  orbit  eccentricity  of  0.3.  The 
configuration  constraint  is  a  fixed  radius  from  the  reference  point.  That  is,  the  swarm  is 
restricted  to  a  surface  that  lies  on  the  sphere  r  =  constant .  This  path  constraint  is  written 

g  =  (fx2+f+r2)  =  r2  (4.25) 

The  goal  is  to  minimize  the  fuel  required  to  meet  the  configuration  constraints 
shown  in  Table  IV-4.  Figure  IV- 13  shows  the  three-dimensional  plot  of  the  solution 
formation.  It  closely  resembles  Figure  III- 1 ,  the  C-W  circular  formation,  but  closer 
inspection  will  reveal  the  subtle  differences.  Figure  IV-14  through  Figure  IV-16  show 
the  projection  of  the  fonnation  in  the  three  orthogonal  planes.  Figure  IV-17  shows  the 
open  loop  controls  required  by  one  of  the  satellites  to  maintain  this  fonnation  for  the 
given  eccentricity.  This  solution  was  found  using  100  nodes  and  fixing  the  final  time  to 
exactly  one  orbit. 
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Cross-Track 


Table  IV-4  Constraints  for  Forced  Circular  Formation 


Constraint 

Normalized 

Lower  and  Upper  Bounds 

States:  rx,ry,rz 

Unconstrained 

States:  rx,ry,rz 

Unconstrained 

States:  m 

[0.1  :  1.0] 

Controls:  T 

[0  :  Unconstrained] 

Time:  x 

%f  =  1.0 

Events:  r(x0)-r(x/) 

[0.0  :  0.0] 

Events:  f(x0)-r(x/) 

[0.0  :  0.0] 

Path:  g  =  rx  +  r2  +  r2 

[1.0  :  1.0] 

Number  of  Nodes 

100 

Reference  Orbit:  e 

0.3 
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Cross-Track  <£'  Radial 


ure  IV- 14  Forced  Circular  Formation  Radial  vs.  Along  Track  Motion 


Radial 


Figure  IV- 15  Forced  Circular  Formation  Cross-Track  vs.  Radial  Motion 


Cross-Track 
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A  comparison  of  Figure  IV- 15  and  Figure  III-3  shows  the  departure  from  the 
classic  C-W  solution.  Figure  IV- 15  apparently  shows  non-planar  motion,  but  itself  is  not 
enough  to  determine  if  the  formation  actually  resides  in  a  plane.  The  answer  lies  in  the 
angular  momentum  vector,  hre/=rxr .  The  direction  of  the  relative  angular  moment 

vector  was  plotted  to  determine  if  the  motion  was  planar.  If  indeed  the  motion  is  planar, 
then  the  direction  of  the  angular  momentum  should  remain  constant.  If  the  formation  is 
not  planar,  then  the  direction  of  the  momentum  vector  will  not  be  constant.  Figure  IV- 18 
shows  this  plot  for  the  classic  C-W  circular  fonnation  described  in  section  III.D.2  and 
demonstrates  that  it  is  constant,  within  numeric  tolerances.  On  the  other  hand,  Figure 
IV- 19  shows  the  same  trace  for  this  fonnation.  It  is  clearly  not  constant  and  therefore 
demonstrates  that  this  formation  is  non-planar.  In  both  figures,  a  line  is  drawn  from  the 
origin  in  the  direction  of  the  unit  vector  associated  with  the  relative  angular  momentum  at 
the  first  time  step.  The  dots  depict  the  tip  of  this  angular  momentum  unit  vector  at  each 
time  step. 


Figure  IV- 18  Relative  Angular  Momentum  Vector  for  the  C-W  Circular  Formation 
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Figure  IV-19  Relative  Angular  Momentum  Vector  for  the  Forced  Circular  Fonnation 


2.  Forced  Projected  Circular  Formation 

This  fonnation  in  another  example  of  a  forced  formation  using  a  well-known 
configuration.  The  configuration  constraint  is  the  same  as  seen  in  section  III.D.3,  except 
for  a  reference  orbit  eccentricity  of  0.3.  The  goal  remains  the  same:  to  minimize  the  fuel 
required  to  meet  the  configuration  constraints  shown  in  Table  IV-5.  Figure  IV-20  shows 
the  three-dimensional  plot  of  the  solution  formation.  Figure  IV-21  through  Figure  IV-23 
show  the  projection  of  this  formation  in  the  three  orthogonal  planes.  Figure  IV-24  shows 
the  open  loop  controls  required  by  one  of  the  satellites  to  maintain  this  formation  for  the 
given  eccentricity.  This  solution  was  also  found  using  100  nodes  and  fixing  the  final 
time  to  exactly  one  orbit. 
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Cross-Track 


Table  IV- 5  Constraints  for  Forced  Projected  Circular  Formation 


Constraint 

Normalized 

Lower  and  Upper  Bounds 

States:  rx,ry,rz 

Unconstrained 

States:  rx,ry,rz 

Unconstrained 

States:  m 

[0.1  :  1.0] 

Controls:  T 

[0  :  Unconstrained] 

Time:  x 

%f  =  1.0 

Events:  r(x0)-r(x/) 

[0.0  :  0.0] 

Events:  f(x0)-r(x/) 

[0.0  :  0.0] 

Path:  g  =  r  +  rz 

[1.0  :  1.0] 

Number  of  Nodes 

100 

Reference  Orbit:  e 

0.3 

1 


Figure  IV-20  Forced  Projected  Circular  Fonnation  for  e  =  0.3 
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Figure  IV -2 1  Forced  Projected  Circular  Formation  Radial  vs.  Along  Track  Motion 


Radial 

Figure  IV-22  Forced  Projected  Circular  Formation  Cross-Track  vs.  Radial  Motion 


E.  VALIDATION  OF  SOLUTIONS 

For  all  the  cases  presented  here,  the  following  steps  were  taken  in  validating  the 
solutions.  The  Lagrangian  was  constructed  in  accordance  with  section  II.B.l  and 
equation  (2.19).  Care  must  be  taken  to  include  the  state  and  control  bounds  in  defining 
the  path  constraints,  g  .  Taking  the  partial  derivatives  of  the  Lagrangian  with  respect  to 
the  controls  yields  the  following  equations: 
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^Tz_ 

Xf  —  T0 

m 

(|)  is  called  the  switching  function  for  the  controls  and  is  governed  by  the  KKT  conditions 
described  in  equation  (2.22).  Figure  IV -25  shows  the  switching  function  for  the  Thrust  in 
each  of  the  three  axes  corresponding  to  the  thrust  profile  shown  in  Figure  IV-24.  Not 
inherently  obvious  is  that  the  switching  function  is  the  same  for  both  circular  and 
projected  circular  fonnations  at  given  values  of  v,  since  the  path  constraints  are  not 
control  dependent. 
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Figure  IV-25  Switching  Functions  for  Control  Thrust 


In  addition  to  the  optimality  analysis  described  above,  the  initial  conditions  from 
the  DIDO  solution  were  numerically  propagated  for  50  orbits.  For  the  natural  solutions, 
the  initial  conditions  were  propagated  with  all  components  of  thrust  set  equal  to  zero.  For 
the  forced  solutions,  the  solution  thrust  profile  was  imported  into  the  propagator.  This 
profile  was  interpolated  to  determine  the  value  for  the  controls  at  each  propagation  time 
step.  The  propagator  was  free  to  determine  its  own  time  steps,  which  did  not  necessarily 
match  the  solution  time  steps.  In  the  process  of  creating  a  propagator,  every  MATLAB 
resident  ODE  solver  was  evaluated.  ODE45  proved  to  have  the  best  combination  of 
accuracy  and  speed.  Additionally,  many  different  interpolation  schemes  were  evaluated. 
A  MATLAB  function  called  POLINT  provided  the  most  accurate  results.  It  was  created 
by  Weideman  and  Reddy  and  implements  the  barycentric  fonnula  from  Henrici’s 
Essentials  of  Numeric  Analysis.  Spline  interpolation  was  the  next  best  in  accuracy,  and 
provided  a  significant  benefit  in  speed  over  POLINT.  Spline  was  used  for  routine 
evaluations  and  POLINT  was  reserved  for  more  detailed  assessments. 
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1.  Natural  Formation  for  e  =  0.5 

Figure  IV-26  and  Table  IV-6  show  that  there  is  no  appreciable  deviation  in 
formation  configuration  after  being  propagated  for  50  orbits. 


Figure  IV-26  Natural  Formation  with  e  =  0.5  Over  50  Orbits 


Table  IV-6  Propagation  Results  for  Natural  Formation  with  e  =  0.5 


State  Variable 

%  Difference  Between 
Initial  and  Final  Values 

8.54  xlO'4 

r 

y 

5.77  xlO'1 

rz 

4.61  xlO'5 

K 

6.67  xlO'1 

1.04x1 0'3 

K 

2.35  xlO'4 
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2.  Natural  Formation  for  e  =  0.3 

Again,  the  initial  conditions  were  numerically  propagated  (with  thrust  equal  to 
zero)  for  50  orbits.  Figure  IV-27  and  Table  IV-7  show  that  there  is  no  significant 
difference  over  this  period. 


Figure  IV-27  Natural  Formation  with  e  =  0.3  Over  50  Orbits 


Table  IV-7  Propagation  Results  for  Natural  Fonnation  with  e  =  0.3 


State  Variable 

%  Difference  Between 
Initial  and  Final  Values 

9.06  xlO'4 

r 

y 

8.71  xlO'2 

rz 

1.97x1  O'4 

K 

5.73  xlO'2 

ry 

1.11  xlO'3 

K 

6.67  xlO'5 
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3.  Natural  Formation  for  e  =  0.7 

Since  solution  three  is  a  reproduction  of  a  previously  known  solution,  the 
validation  of  the  fonnation  configuration  can  be  found  in  Reference  17. 

4.  Forced  Circular  Formation 

The  validation  of  this  result  included  propagating  the  initial  conditions 
numerically  for  50  orbits,  now  subject  to  the  thrust  profile  shown  in  Figure  IV-17.  Figure 
IV-28  shows  the  three-dimensional  motion  of  the  formation  over  this  time  period.  Table 
IV-8  also  demonstrates  that,  as  desired,  there  is  no  appreciable  digression  in  the 
formation  configuration  over  50  orbits. 
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Figure  IV-28  Forced  Circular  Formation  over  50  Orbits 
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Table  IV- 8  Propagation  Results  for  Forced  Circular  Formation 


State  Variable 

%  Difference  Between 
Initial  and  Final  Values 

C 

0.23 

ry 

1.51 

C 

0.04 

K 

0.22 

0.80 

K 

0.11 

One  unique  consequence  of  imposing  the  circular  fonnation  constraint  is  that  it 
provides  another  way  to  verify  the  results.  The  forced  circular  is  defined  by  r  =  constant 

which  means  r  =  0 .  This  implies  r 2  =  r  ■  r  =  const ,  which  gives  (r2 )  =  2(r •  r)  =  0 , 

followed  naturally  by 

— (r  r)  =  r  r  +  r  r  =  0  (4.32) 

dt 

Substituting  the  EOM  for  r  above  provides  an  independent  check  that  the  solution  is 
indeed  a  circular  formation. 

5.  Forced  Projected  Circular  Formation 

As  with  the  previous  solutions,  the  initial  conditions  were  propagated  numerically 
for  50  orbits,  subject  to  the  thrust  profile  shown  in  Figure  IV-24.  Figure  IV-29  shows  the 
three-dimensional  motion  of  the  formation  over  this  time.  Table  IV-9  also  details  the 
errors  in  the  fonnation  configuration  after  50  orbits. 
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Cross-Track 


Figure  IV-29  Forced  Projected  Circular  Formation  over  50  orbits 


Table  IV-9  Propagation  Results  for  Forced  Projected  Circular  Formation 


State  Variable 

%  Difference  Between 
Initial  and  Final  Values 

C 

0.06 

ry 

1.33 

rz 

0.01 

K 

0.11 

c 

0.47 

K 

0.02 

61 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


62 


V.  UNRESOLVED  ISSUES 


A.  OPTIMAL  PERIOD  VS.  ORBITAL  PERIOD 

For  Natural  or  Zero-Thrust  solutions,  the  optimal  period,  (x/-x0J,  is  exactly 

equal  to  the  orbital  period.  However,  once  thrust  is  required  to  maintain  the  formation 
the  optimal  period  may  or  may  not  equal  the  orbital  period.  One  way  to  resolve  this  is  to 
force  the  solution  period  to  be  equal  to  the  orbital  period.  However,  this  fix  raises  several 
questions. 

Table  V-l  shows  the  cost  associated  with  two  different  solutions  to  the  Forced 
Circular  Formation  with  e  =  0.3  and  assumes  a  100  kg  satellite,  with  1  =  1000  sec  ,  and 

the  TU  for  a  perigee  altitude  of  1000  km.  Solution  1A  is  detailed  in  section  IV.D.l, 
while  Solution  2A  is  not  shown  in  the  prior  chapters.  Table  V-2  shows  the  cost 
associated  with  two  different  solutions  to  the  Forced  Projected  Circular  Fonnation  using 
the  same  assumptions  described  above.  Solution  IB  is  represented  in  section  IV.D.2, 
while  Solution  2B  is  not  shown. 


Table  V-l  Differing  Costs  for  Forced  Circular  Solutions 


Forced  Circular  Solution  1A 

Cost:  t/^Tj  —  r0) 

1.98  xlO'7  Newtons/sec 

Cost:  %  Mass 

2.34  xlO'5  % 

Optimal  Period:  xf 

1 .0  orbital  periods 

Forced  Circular  Solution  2A 

O 

o 

VI 

r-t- 

1 

1.80  xlO'7  Newtons/sec 

Cost:  %  Mass 

1.00  xlO'5  % 

Optimal  Period:  x/ 

1.0005  orbital  periods 
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Table  V-2  Differing  Costs  for  Forced  Projected  Circular  Solutions 


Forced  Projected  Circular  Solution  IB 

Cost:  T/(Tf  -t0) 

2.87  xlO'7  Newtons/sec 

Cost:  %  Mass 

3.39  xlO'5  % 

Optimal  Period:  xf 

1 .0  orbital  periods 

Forced  Projected  Circular  Solution  2B 

Cost:  T/[Tf  -  r0) 

2.85  xlO'7  Newtons/sec 

Cost:  %  Mass 

3.38  xlO'5  % 

Optimal  Period:  xf 

1 .002  orbital  periods 

In  both  cases,  the  second  solution  has  the  lower  total  cost.  However,  the  optimal 
period  of  these  solutions  is  not  equal  to  the  orbital  period.  The  precise  explanation  is 
unknown.  These  solutions  may  be  periodic  as  a  fonnation,  but  not  periodic  with  respect 
to  the  reference  point  since  the  period  of  the  reference  point  (or  satellite)  is,  by  definition, 
equal  to  the  orbital  period.  This  mismatch  in  periods  presents  difficulty  in  visualizing  or 
plotting  the  configuration  or  relative  motion.  If  ignored,  the  formation  appears  to  drift 
away  from  the  reference  point  and  therefore  not  stay  together.  One  way  to  overcome  this 
visualization  problem  is  to  display  the  position  of  one  swann  satellite  against  another 
swarm  satellite.  Over  the  course  of  the  optimal  period,  the  distance  should  remain 
bounded  and  repeat  itself. 

B.  SWITCHING  FUNCTION  :  DERIVED  VS.  DIDO  SOLUTION 

Normally,  (see  equations  (4.26)  through  (4.31)  )  is  equal  to  zero  for  an 
dT  ' 

optimal  solution.  This  allows  the  values  for  the  switching  function,  ()> ,  to  be  derived  and 
calculated.  At  the  same  time,  DIDO  calculates  the  values  for  the  switching  function 
directly.  Figure  IV-25  shows  the  switching  function  as  calculated  by  DIDO,  which 

agrees  with  the  expected  results.  The  concern  is  that  the  value  for  is  not  0.0  but 

dT 
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exactly  0.5.  The  resulting  switching  functions  have  the  identical  shape  as  those  in  Figure 
IV-25,  but  are  offset  by  negative  0.5  in  the  vertical  axis.  This  discrepancy  in  the 
switching  function  requires  further  analysis. 

C.  J2  PERTURBATIONS 

Work  was  started  on  implementing  the  effects  of  earth  oblateness  or  J2  in  the 
relative  frame,  but  due  to  time  constraints  was  not  completed. 

1.  Linear  J2  terms. 

The  following  Linear  equations  for  J2  effects  in  the  relative  frame  came  from  an 
unpublished  paper  by  I.M.  Ross.  Using  the  same  coordinate  system  as  Figure  II- 1  and 
assuming  a  circular  reference  orbit,  the  J2  effects  can  be  represented  as  relative 
accelerations. 


12 sin"  z'sin2  05 

-4  sin2  z  sin  205 

-4  sin  2z  sin  05 

rx 

aj2  n  Jr 

-4 sin2  z'sin 205 

1  +  sin2  z  (2 -7 sin2  05) 

sin  2z  cos  05 

ry 

-4  sin  2  z  sin  05 

sin  2  z  cos  05 

3-sin2  z'(2  +  5sin2  05) 

_L_ 

(5.1) 
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and  05  =  03  + v 
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Non-linear  J2  terms. 


20 

By  comparison,  the  non-linear  J2  terms-  in  the  inertial  reference  frame  are 
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where 


JD 


-3J1\xR2e 

2 


(5.6) 


In  order  to  get  the  relative  accelerations  resulting  from  J2  into  the  relative  frame,  several 
transfonnations  must  be  defined.  To  transform  perifocal  to  Earth-centered  inertial  (ECI): 


cQcco-sQscocz 

—cQs(0—sQc(Oci 

sClsi 

IJK 

— 

sQc(0  +  cQ.s(x)ci 

sQsm+cQccocz 

- cQsi 

PQW 

sCdsi 

cQdsi 

ci 

with  c  and  s  representing  sin  and  cosine  of  the  appropriate  angles.  The  transpose  of 
above  is  used  to  convert  from  ECI  to  perifocal.  Transforming  the  formation  reference 
frame  to  perifocal  is  done  by 


COSY 

-sinv 

0" 

“  PQW~ 

0 

RSW  _ 

— 

sinv 

COSY 

0 

0 

1 

with  the  transpose  used  for  moving  in  the  other  direction. 


The  algorithm  used  to  calculate  J2  from  the  relative  positions  of  the  swarm 
satellites  begins  with  converting  relative  positions  into  inertial  positions. 


~  PQW  ~ 

R,-e/  + 

PQW 

"  RSW " 

IJK 

IJK 

PQW 

(5.9) 


The  next  step  is  to  calculate  the  accelerations  in  the  inertial  frame  according  to  equations 
(5.3)  through  (5.5).  These  inertial  accelerations  can  now  be  transformed,  using  equation 
(4.15)  into  the  relative  accelerations.  Of  note  is  the  use  of  a  static  reference  orbit.  That 
is,  the  reference  orbit  is  assumed  to  be  unperturbed. 

3.  Comparison  Between  Sets  of  Equations. 

Initial  results  from  a  comparison  between  the  linear  equations  and  the  non-linear 

equations  were  encouraging.  The  comparison  was  done  for  a  circular  formation,  with 

e  =  0  at  three  different  inclinations:  28.5°,  45°,  and  63.4°.  Figure  V-l  to  Figure  V-3  show 

the  perturbation  accelerations  due  to  J2  effects  for  all  three  inclinations,  nonnalized  by  the 
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scaling  scheme  detailed  in  section  II. B. 2.  Figure  V-4  to  Figure  V-6  show  the  difference 
between  equation  (5.1)  and  equations  (5.3)  through  (5.5)  for  each  inclination.  Note  that 
the  differences  are  not  random,  but  seem  to  follow  some  pattern.  This  is  the  effect  of  the 
higher  order  tenns  not  present  in  the  linear  equations. 


Figure  V-l  Relative  Accelerations  Due  to  J2  at  28.5°  Inclination 
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Relative  J_  Accelerations  Relative  J0  Accelerations 
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Figure  V-2  Relative  Accelerations  Due  to  J2  at  45°  Inclination 
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Figure  V-3  Relative  Accelerations  Due  to  J2  at  63.4°  Inclination 
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VI.  CONCLUSIONS  AND  FUTURE  WORK 


A.  CONCLUSIONS 

This  work  for  this  thesis  started  as  a  question:  Can  optimal  periodic  control 
(OPC)  theory  be  applied  to  satellite  formation  design  and  control?  Using  the 
fundamentals  of  optimal  control  theory,  a  framework  was  developed  that  captures  the 
essence  of  designing  and  controlling  spacecraft  formations.  This  framework  is  used  to 
articulate  a  variety  of  formations  including  the  notion  of  an  aperiodic  formation.  Based 
on  a  deliberate  formulation,  including  mass  as  a  state  variable,  it  was  shown  that  the 
numerical  approach  can  easily  handle  nonlinearities.  Additionally,  fonnations  were 
presented  that  require  active  control  along  with  their  corresponding  thrust  profile. 

This  thesis  lends  credence  to  the  notion  of  numerically  searching  for  minimum- 
fuel  formation  configurations  for  spacecraft  swarm  subject  to  arbitrary  nonlinear 
dynamics.  Thus,  practical  formations  may  be  designed  and  controlled  using  this  method. 
The  foundation  for  applying  OPC  to  satellite  formations  has  been  completed.  Future 
work  can  build  on  this  foundation  in  increasing  the  complexity  of  the  models.  A 
substantial  amount  of  work  needs  to  be  done  in  developing  a  more  complete  model  of 
satellite  motion  and  perturbations.  The  framework  developed  here  can  be  readily  used 
with  more  robust  dynamical  models  and  should  produce  very  interesting  results. 

B.  OPPORTUNITIES  FOR  FUTURE  WORK 

In  addition  to  the  topics  identified  in  the  previous  chapters,  opportunities  for 
future  work  are  available  in  the  following  areas. 

In  the  course  of  implementation  for  Chapter  IV,  an  additional  “state”  was  added 
to  simplify  the  equations.  This  state  variable  was  v  with  its  corresponding  derivative,  v  . 
This  required  adding  dynamical  constraints  for  this  state.  One  opportunity  for  future 
work  would  be  to  include  v  in  the  path  constraints  and  rewrite  the  dynamical  equations 
to  remove  explicit  dependence  on  v .  This  may  possibly  increase  numerical  accuracy 
when  implementing  the  model  in  DIDO. 
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Another  possibility  for  follow  on  research  is  an  examination  of  the  Nd)B  term  .  It 
may  be  possible  to  continue  working  in  the  formation  reference  frame,  while  adding  the 

J2  effects  to  this  term  by  including  —  anci  <10)- 

dt  dt  dt 

A  departure  from  the  relative  reference  frame  would  be  an  area  for  future  work. 
The  creation  of  a  model  in  an  inertial  frame  and  the  use  of  the  full  non-linear  equations 
will  create  a  high  fidelity  model,  and  will  allow  the  addition  of  J2  effects  and  any  other 
Earth-centered  effects  to  be  included  more  readily.  Translating  the  formation  constraints 
from  the  relative  frame,  where  they  are  intuitive,  to  the  inertial  frame  should  be 
straightforward  but  tedious.  If  the  model  is  defined  in  inertial  space,  the  deployment  and 
reconfiguration  problems  becomes  a  simple  change  from  one  orbit  to  another  orbit.  The 
difficulty  lies  in  specifying  the  constraints,  if  any,  on  the  reconfiguration  relative  motion 
and  visualization. 


Another  consideration,  offered  by  Dr.  Terry  Alfriend,  is  to  write  the  path 
constraint  as  follows,  and  minimize  total  thrust. 


f  V 


\  R«r  J 


f  V 

ry 

v  R,  ef  J 


+ 


R 


ref 


=  constant 


(6.1) 


If  the  reference  orbit  is  circular,  then  Rref  is  constant  and  becomes  a  simple  scaling 

factor,  which  was  completed  previously.  The  above  constraint  was  not  addressed  for  an 
elliptical  reference  orbit,  where  Rref  is  no  longer  a  constant. 

Other  unexplored  areas  include  the  possibility  of  near-circular  relative  formations. 
If  a  small  deviation  is  allowed  in  the  path,  it  may  lower  the  cost  of  the  formation.  In 
addition,  specific  missions  require  specific  configuration  constraints.  One  very 
interesting  area  of  research  would  be  to  identify  and  catalog  various  configuration 
constraints  according  to  mission. 
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